##############################
##############################
##### Lillios, Moskowitz, Resch
##### Prejudice and Turnout
##### Gov 2001
##### 5 May 2014
##############################
##############################

# clear workspace
rm(list=ls())

# set working directory
setwd ("#insert wd link here")

# open necessary libraries
library(foreign)
library(psych)
library(stargazer)
library(Matching)
library(ebal)
library(Zelig)
library(survey)
library(sandwich)


# read in the data set
anes <- read.dta("NES Time Series 2008.dta",
                 convert.factors = F)

##########################
###### DATA CLEANING #####
##########################

##### functions used for data cleaning ####
### recode from 0 to 1 ###
### this function recodes variables on a 0-1 scale as is done in 
### Krupnikov and Piston (2014)
recoder <- function(old.var){
  min.old <- min(old.var, na.rm=T)
  max.old <- max(old.var, na.rm=T)
  new.var <- (old.var - min.old) / (max.old - min.old)
  return(new.var)
}


##### dependent variable #####
### turnout variable
anes$turnout <- anes$V085036x
anes$turnout[anes$turnout == -2] <- NA
prop.table(table(anes$turnout))


##### explanatory variables #####
### party id
anes$party.id <- anes$V083097
# refused recoded as NA
anes$party.id[anes$party.id==-9] <- NA
# don't know (=-8), other (=4), no pref (=5) recoded as independent (=3)
anes$party.id[anes$party.id %in% c(-8,4,5,3)] <- 3
# republicans recoded =5
anes$party.id[anes$party.id==2] <- 5
# if indep, no pref, other, dk and closer to dems recode party id =2
anes$party.id[anes$party.id==3 & anes$V083098b==5] <- 2 
# if indep, no pref, other, dk and closer to reps recode party id =4# most
anes$party.id[anes$party.id==3 & anes$V083098b==1] <- 4 
# if dem and identify as strong dem, recode party id =0
anes$party.id[anes$party.id==1 & anes$V083098a==1] <- 0 
# if rep and identify as strong rep, recode party id =6
anes$party.id[anes$party.id==5 & anes$V083098a==1] <- 6 
# create a party id standardized from 0 to 1
anes$party.id.0.1 <- recoder(anes$party.id)
# validate
table(anes$party.id)
table(anes$party.id.0.1)

### party id dummy variables
anes$strong.dem[anes$party.id==NA] <- NA
anes$weak.dem[anes$party.id==NA] <- NA
anes$lean.dem[anes$party.id==NA] <- NA
anes$ind[anes$party.id==NA] <- NA
anes$lean.rep[anes$party.id==NA] <- NA
anes$weak.rep[anes$party.id==NA] <- NA
anes$strong.rep[anes$party.id==NA] <- NA
# strong dem
anes$strong.dem[!is.na(anes$party.id)] <- 0
anes$strong.dem[anes$party.id==0] <- 1
# weak dem
anes$weak.dem[!is.na(anes$party.id)] <- 0
anes$weak.dem[anes$party.id==1] <- 1
# lean dem
anes$lean.dem[!is.na(anes$party.id)] <- 0
anes$lean.dem[anes$party.id==2] <- 1
# independent
anes$ind[!is.na(anes$party.id)] <- 0
anes$ind[anes$party.id==3] <- 1
# lean republican
anes$lean.rep[!is.na(anes$party.id)] <- 0
anes$lean.rep[anes$party.id==4] <- 1
# weak republican
anes$weak.rep[!is.na(anes$party.id)] <- 0
anes$weak.rep[anes$party.id==5] <- 1
# strong republican
anes$strong.rep[!is.na(anes$party.id)] <- 0
anes$strong.rep[anes$party.id==6] <- 1

### partisan strength
# anes$party.id <- anes$V083097
# anes$strong.partisan <- NA
# anes$strong.partisan[anes$V083098a == 1] <- 1
# anes$strong.partisan[anes$V083098a %in% c(5,-1)] <- 0 
# # validate our results
# prop.table(table(anes$strong.partisan))


### prejudice scale
### we construct a difference measure for both the lazy" question and the 
### "unintelligent" question by subtracting the score the respondent gave 
### Blacks from the score for Whites... The number 6 is added to all points on 
### the scale, resulting in a scale that ranges from 0 to 12. All points on 
### the scale are then divided by 12. (see Krupnikov and Piston 2014, p. 5).
## DM: While they claim to subtract the score for Blacks from the score for 
## Whites, it appears that they actually subtract the score for Whites from
## the score for Blacks.
# assign NA to all values < 1 for variables used in the prejudice scale.
prej.varlist <- colnames(with(anes,cbind(V083207a,V083207b,V083208a,V083208b)))
for(i in 1:length(prej.varlist)){
  anes[!is.na(anes[,prej.varlist[i]]) < 1,prej.varlist[i]] <- NA
}
# create the prejudice scale
lazy.scale <- with(anes,(V083207b - V083207a + 6)/12)
describe(lazy.scale)
intelligent.scale <- with(anes, (V083208b - V083208a + 6)/12)
describe(intelligent.scale)
anes$prejudice <- (lazy.scale + intelligent.scale) / 2
describe(anes$prejudice)

# create a p prejudice scale from the post-election data
prej.post.varlist <- colnames(with(anes,cbind(V085174a,V085174b,V085175a,V085175b)))
for(i in 1:length(prej.varlist)){
  anes[!is.na(anes[,prej.varlist[i]]) < 1,prej.varlist[i]] <- NA
}
lazy.scale.post <- with(anes,(V085174b - V085174a + 6)/12)
describe(lazy.scale)
intelligent.scale <- with(anes, (V083208b - V083208a + 6)/12)
describe(intelligent.scale)
anes$prejudice <- (lazy.scale + intelligent.scale) / 2
describe(anes$prejudice)

### age
# in years
anes$age <- anes$V083215x
anes$age[anes$age < 18] <- NA
describe(anes$age)
# recoded 0-1
anes$age.0.1 <- recoder(anes$age)
describe(anes$age.0.1)

### education
# highest grade
anes$educ <- anes$V083217
anes$educ[anes$educ < 0] <- NA
describe(anes$educ)
# recoded 0-1
anes$educ.0.1 <- recoder(anes$educ)
describe(anes$educ.0.1)

### gender: female coded =1
anes$female <- anes$V081101 - 1
table(anes$female)

### income
anes$income <- anes$V083249
anes$income[anes$income < 0] <- NA
anes$income.0.1 <- recoder(anes$income)
describe(anes$income.0.1)

### currently married
anes$married <- NA
anes$married[anes$V083216x %in% 2:6] <- 0
anes$married[anes$V083216x==1] <- 1
table(anes$married)

### south
states.in.south <- c("FL","SC","AL","MS","GA","LA","TX","VA","AR","TN","NC")
anes$south <- 0
anes$south[anes$V081201a %in% states.in.south] <- 1
table(anes$south)

### interest
# old version of the political interest variable:
old.interest <- NA
# most of the time coded =1
old.interest[anes$V085072==1] <- 1
# some of the time coded =.67
old.interest[anes$V085072==2] <- 0.67
# only now and then coded =.33
old.interest[anes$V085072==3] <- 0.33
# hardly at all coded  =0
old.interest[anes$V085072==4] <- 0
# new version of the political interest variable:
new.interest <- NA
# extremely closely coded =1
new.interest[anes$V085073a==1] <- 1
# very closely coded = 1
new.interest[anes$V085073a==2] <- 1
# moderately closely coded =.67
new.interest[anes$V085073a==3] <- .67
# slightly closely coded =.33
new.interest[anes$V085073a==4] <- .33
# not closely at all coded =0
new.interest[anes$V085073a==5] <- 0
# create a combined interest variable from both versions
anes$interest <- NA
anes$interest <- old.interest
anes$interest[is.na(old.interest)] <- new.interest[is.na(old.interest)]
table(anes$interest)

### contacted by party
# DM: possible condition: & anes$V085025a %in% c(1,5,6)
anes$contact <- NA
anes$contact[anes$V085025==1] <- 1
anes$contact[anes$V085025==5] <- 0
table(anes$contact)

### length of resident in R's community
# DM: THEY DO NOT CODE 77 PROPERLY FOR V083266a
# 77 INDICATES THEY HAVE LIVED THERE ENTIRE LIFE NOT 77 YEARS
# in years 
residence <- NA
residence[anes$V083266a>=0 & !is.na(anes$V083266a)] <- 
  anes$V083266a[anes$V083266a >= 0 & !is.na(anes$V083266a)]
# recoded 0 to 1
anes$residence.0.1 <- recoder(residence)
# validate
describe(anes$residence.0.1)

### external efficacy
## public officials do (not) care
# version C (V083079c): public officials don't care about what people think
officials.care <- NA
# agree strongly coded =0
officials.care[anes$V083079c==1] <- 0
# agree somewhat coded =.25
officials.care[anes$V083079c==2] <- .25
# neither agree nor disagree coded =.50
officials.care[anes$V083079c==3] <- .50
# disagree somewhat coded =.75
officials.care[anes$V083079c==4] <- .75
# disagree strongly coded =1
officials.care[anes$V083079c==5] <- 1
# version D (V083080c): officials care what people like you think
# a great deal coded =1
officials.care[anes$V083080c==1] <- 1
# a lot coded =.75
officials.care[anes$V083080c==2] <- .75
# a moderate amount coded =.50
officials.care[anes$V083080c==3] <- .50
# a little coded =.25
officials.care[anes$V083080c==4] <- .25
# not at all coded =0
officials.care[anes$V083080c==5] <- 0
# validate
table(officials.care)
## people have a (no) say in what government does
# version C (V083079d): People don't have any say about what government does.
say.in.gov <- NA
# agree strongly coded =0
say.in.gov[anes$V083079d==1] <- 0
# agree somewhat coded =.25
say.in.gov[anes$V083079d==2] <- .25
# neither agree nor disagree coded =.50
say.in.gov[anes$V083079d==3] <- .50
# disagree somewhat coded =.75
say.in.gov[anes$V083079d==4] <- .75
# disagree strongly coded =1
say.in.gov[anes$V083079d==5] <- 1
# version D (V083080d): How much can people affect what government does? 
# a great deal coded =1
say.in.gov[anes$V083080d==1] <- 1
# a lot coded =.75
say.in.gov[anes$V083080d==2] <- .75
# a moderate amount coded =.50
say.in.gov[anes$V083080d==3] <- .50
# a little coded =.25
say.in.gov[anes$V083080d==4] <- .25
# not at all coded =0
say.in.gov[anes$V083080d==5] <- 0
# validate
table(say.in.gov)
## combine officials.care and say.in.gov for external efficacy index
anes$external.efficacy <- (officials.care + say.in.gov) / 2
table(anes$external.efficacy)

### does R think presidential race will be close
anes$close.race <- NA
anes$close.race[anes$V083074==1] <- 1
anes$close.race[anes$V083074==5] <- 0
table(anes$close.race)

### affect for candidate
obama.therm <- anes$V083037a
obama.therm[obama.therm < 0] <- NA
mccain.therm <- anes$V083037b
mccain.therm[mccain.therm < 0] <- NA
anes$affect <- obama.therm - mccain.therm
describe(anes$comb.interest[!is.na(anes$turnout)])

##### VARIABLES FOR SUBSETTING AND WEIGHTS #####
### race
# white
anes$white <- NA
anes$white[!is.na(anes$V083251a) & anes$V083251a>0] <- 0
anes$white[anes$V083251a==50] <- 1
table(anes$white)
# non-hispanic white
anes$nh.white <- anes$white
anes$nh.white[anes$white==1 & (anes$V083255==1 | anes$V083256==1)] <- 0
table(anes$nh.white)

# create a variable for strong dem and nh white
anes$strong.dem.nh.white <- NA
anes$strong.dem.nh.white[!is.na(anes$strong.dem) & !is.na(anes$nh.white)] <- 0
anes$strong.dem.nh.white[anes$strong.dem==1 & anes$nh.white==1] <- 1

### survey weights
anes$weight <- anes$V080102

### psu
anes$psu <- anes$V081205

####################
##### ANALYSIS #####
####################

### MODELS 1-6 in TABLE 1 ###
# variables for models 1,2,3,4
model.1234.X <- as.matrix(with(anes,cbind(
  turnout,prejudice,age.0.1,educ.0.1,female,income.0.1,married,residence.0.1,
  south,external.efficacy,interest,contact,close.race,weight,psu)))
colnames(model.1234.X) <- c("Turnout","Prejudice","Age","Education","Female",
                            "Income","Married","Years in residence","South",
                            "Efficacy","Interest","Contacted by party",
                            "Perceive election to be close","weight","psu")
# variables for model 5
model.5.X1 <- as.matrix(with(anes,cbind(
  turnout,prejudice,strong.dem,prejudice*strong.dem,age.0.1,educ.0.1,female,
  income.0.1,married,residence.0.1,south,external.efficacy,interest,
  contact,close.race,weight,psu)))
colnames(model.5.X1) <- c("Turnout","Prejudice","Strong Dem",
                         "Strong Dem x Prej","Age","Education","Female",
                         "Income","Married","Years in residence","South",
                         "Efficacy","Interest","Contacted by party",
                         "Perceive election to be close","weight","psu")
# variables for model 6
model.6.X1 <- as.matrix(with(anes,cbind(
  turnout,prejudice,strong.dem,prejudice*strong.dem,weight,psu)))
colnames(model.6.X1) <- c("Turnout","Prejudice","Strong Dem",
                         "Strong Dem x Prej","weight","psu")

# model 1: subset to strong dems (and whites)
model.1.X <- na.omit(subset(model.1234.X,anes$strong.dem==1 & anes$nh.white==1))
# model 2: subset to weak/lean dems (and whites)
model.2.X <- na.omit(subset(model.1234.X,(anes$weak.dem==1 | anes$lean.dem==1) &
                      anes$nh.white==1))
# model 3: subset to independents (and whites)
model.3.X <- na.omit(subset(model.1234.X,anes$ind==1 & anes$nh.white==1))
# model 4: subset to republicans (and whites)
model.4.X <- na.omit(subset(model.1234.X,(anes$strong.rep==1 | anes$weak.rep==1 |
                                    anes$lean.rep==1) & anes$nh.white==1))
# model 5: all whites
model.5.X <- na.omit(subset(model.5.X1,anes$nh.white==1))
# model 6: all whites
model.6.X <- na.omit(subset(model.6.X1,anes$nh.white==1))

# create a function for the matrices above to run logit
# logit.Krup.Pist <- function(data.matrix){
#   glm(data.matrix[,1] ~ data.matrix[,2:(NCOL(data.matrix)-1)], 
#       family = "binomial"(link = "logit"),
#       weight=data.matrix[,NCOL(data.matrix)])
# }
logit.Krup.Pist <- function(data.matrix){
  glm(data.matrix[,1] ~ data.matrix[,2:(NCOL(data.matrix)-2)], 
      family = "binomial"(link = "logit"),
      weight=data.matrix[,NCOL(data.matrix)-1])
}
# run the function for each model
results.mod1 <- logit.Krup.Pist(model.1.X)
results.mod2 <- logit.Krup.Pist(model.2.X)
results.mod3 <- logit.Krup.Pist(model.3.X)
results.mod4 <- logit.Krup.Pist(model.4.X)
results.mod5 <- logit.Krup.Pist(model.5.X)
results.mod6 <- logit.Krup.Pist(model.6.X)

# weave standard errors to replicate K-P's from stata
weave.se <- function(model){
  sqrt(diag(weave(model)))
}
results.list <- list(results.mod1,results.mod2,results.mod3,results.mod4,results.mod5,results.mod6)
se.list <- lapply(results.list, weave.se)

# output clean results to LaTeX
cov.labels <- c("Prejudice","Strong Dem",
                "Strong Dem x Prej","Age","Education","Female",
                "Income","Married","Years in residence","South",
                "Efficacy","Interest","Contacted by party",
                "Perceive election to be close")
stargazer(results.list,se=se.list,
          title=c("Prejudice and White turnout in the 2008 election, by partisan group (ANES data)"),
          covariate.labels=cov.labels,dep.var.labels=c("Turnout"),column.sep.width="2pt")

###############################
###############################
### quantities of interest ####
###############################
###############################

############################################
### turnout by prejudice for strong dems ###
############################################

prej.vals <- as.numeric(names(table(model.1.X[,"Prejudice"])))
prej.vals.freq <- table(model.1.X[,"Prejudice"])
prej.vals.weightfreq <- tapply(model.1.X[,"weight"],model.1.X[,"Prejudice"],sum)
turnout.by.prej <- prop.table(table(model.1.X[,"Prejudice"],model.1.X[,"Turnout"]),1)
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/turnout_by_prejudice.pdf",
    width=7, height=5)
par(mar=c(4.1,4.1,1.1,2.1))
plot(NA, xlim=c(0,1), ylim=c(0,1.05),
     #main="Turnout by Prejudice for White, Strong Dems",
     ylab="Turnout", 
     xlab="Prejudice", axes=F)
points(x=prej.vals, y=turnout.by.prej[,2], pch=19, type="b")
axis(2, labels=T, lwd.ticks=T)
#text(x=prej.vals,y=turnout.by.prej[,2], labels=round(turnout.by.prej[,2],2), cex=.8, pos=1)# adj=c(-.5,-.5),)
text(x=prej.vals[c(1:3,5,9,12,13)],y=turnout.by.prej[c(1:3,5,9,12,13),2], labels=round(turnout.by.prej[c(1:3,5,9,13,13),2],2), cex=.8, pos=1)
text(x=prej.vals[c(4,6,10)],y=turnout.by.prej[c(4,6,10),2], labels=round(turnout.by.prej[c(4,6,10),2],2), cex=.8, pos=3)
text(x=prej.vals[c(7,11)],y=turnout.by.prej[c(7,11),2], labels=round(turnout.by.prej[c(7,11),2],2), cex=.8, pos=4)
text(x=prej.vals[8],y=turnout.by.prej[8,2], labels=round(turnout.by.prej[8,2],2), cex=.8, pos=2)
par(new = TRUE)
plot(0, axes = F, type="l", bty="n", xlab="", ylab="", xlim=c(0,1), 
     ylim=c(0,1.5*max(prej.vals.freq)))
segments(x0=prej.vals, y0=0, y1=prej.vals.freq, col="gray48", lwd=4)
text(x=prej.vals, y=prej.vals.freq, labels=prej.vals.freq, pos=3, cex=.8, col="gray48")
axis(1, labels=F, lwd.ticks=0)
#mtext(round(prej.vals,2), side=1, at=prej.vals, cex=.6)
mtext(seq(0,1,.25), side=1, at=seq(0,1,.25), cex=.8)
graphics.off()


###############
### model 1 ###
###############
### distribution of prejudice ###
# create a histogram of weighted and unweighted frequencies of prejudice
distr.prej.plotter <- function(data,subtitle=""){
  # calculate the frequencies
  temp.data <- as.data.frame(data)
  x.vals <- as.numeric(names(table(temp.data$Prejudice)))
  x.vals.freq <- table(temp.data$Prejudice)
  sum.weights <- tapply(temp.data$weight,temp.data$Prejudice,sum)
  # plot the unweighted and weighted frequencies
  plot(0, axes = F, 
       type="l", bty="n", xlim=c(0,1), 
       ylim=c(0,max(c(x.vals.freq,sum.weights))),
       #main=paste(subtitle,": Distribution of Prejudice"),
       xlab="", ylab="Frequency")
  segments(x0=x.vals-.006, y0=0, y1=x.vals.freq, col="gray48", lwd=5)
  axis(2)
  par(new = TRUE)
  plot(0, axes = F, type="l", bty="n", xlab="", ylab="", xlim=c(0,1), 
       ylim=c(0,max(c(x.vals.freq,sum.weights))))
  segments(x0=x.vals+.006, y0=0, y1=sum.weights, col="black", lwd=5)
  axis(1, labels=F, lwd.ticks=0)
  #mtext(round(x.vals,2), side=1, at=x.vals, cex=.6)
  mtext(seq(0,1,.25), side=1, at=seq(0,1,.25), cex=.8)
  mtext("Prejudice", side=1, line=1)
  mtext("Unweighted frequencies in gray", side=1, line=2, col="gray48", cex=.8)
  mtext("Weighted frequencies in black", side=1, line=3, col="black", cex=.8)
}
# run the function
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/mod1_distr_prejudice.pdf",
    width=7, height=4)
par(mar=c(4.1,4.1,1.1,2.1))
distr.prej.plotter(data=model.1.X,subtitle="Model 1")
graphics.off()

### predicted probabilities ###
# create a function to plot predicted probabilities
pp.plotter.logit <- function(model,x,interaction=NA,Xc,subtitle="",seed=02143,n=100000){
  # require necessary libraries
  require(mvtnorm)
  # get necessary objects from our glm object
  coefs <- model$coefficients
  #   varcov <- vcov(model)
  varcov <- weave(model)
  # get our simulated betas (beta tildes)
  set.seed(seed)
  sim.betas <- rmvnorm(n = n,
                       mean = coefs,
                       sigma = varcov)
  temp.data <- as.data.frame(model$model[,2])
  x.vals <- as.numeric(names(table(temp.data$Prejudice)))
  x.vals.freq <- table(temp.data$Prejudice)
  # expected values by variable of interest
  pred.probs.mat <- matrix(data=NA, ncol=length(x.vals), nrow=n)
  # for non-interaction models
  if(is.na(interaction)){
    for(j in 1:length(x.vals)){
      new.Xc <- Xc
      new.Xc[x] <- x.vals[j]
      pred.probs.mat[,j] <- (1 / (1 + exp(-new.Xc %*% t(sim.betas))))
    }
  }
  # for interaction models
  if(!is.na(interaction)){
    for(j in 1:length(x.vals)){
      new.Xc <- Xc
      new.Xc[x] <- x.vals[j]
      new.Xc[interaction] <- x.vals[j]
      pred.probs.mat[,j] <- (1 / (1 + exp(-new.Xc %*% t(sim.betas))))
    }
  }
  # predicted probabilities by variable of interest
  pred.prob <- apply(pred.probs.mat,2,mean)
  pred.prob.ci.lb <- apply(pred.probs.mat,2,function(z){quantile(z,.025)})
  pred.prob.ci.ub <- apply(pred.probs.mat,2,function(z){quantile(z,.9725)})
  # plot the predicted probabilities
  plot(x=x.vals,y=pred.prob,
       #main=paste(subtitle,": Predicted Probability of Turnout for Strong Democrats by",x),
       cex.main=0.7,ylab="Predicted Probability of Turnout", 
       ylim=c(0,1),
       xlim=c(0,1),
       xlab=x, cex.lab=.8, cex.axis=.8, pch=19, cex=.8, xaxp=c(0,1,4))
  segments(x0=x.vals, y0=pred.prob.ci.lb,y1=pred.prob.ci.ub)
  # plot the unweighted frequencies
  par(new = TRUE)
  plot(0, axes = F, type="l", bty="n", xlab="", ylab="", xlim=c(0,1), 
       ylim=c(0,2*max(x.vals.freq)))
  segments(x0=x.vals, y0=0, y1=x.vals.freq, col="gray48", lwd=5)
  text(x=x.vals, y=x.vals.freq, labels=x.vals.freq, pos=3, cex=.7, col="gray48")
}

# use the function to plot the predicted probabilities of turnout for values of prejudice
# set all variables to their mean
Xc1 <- c(1,apply(model.1.X[,2:13],2,mean))
# plot the predicted probability of prejudice
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/mod1_pred_probs.pdf",
    width=7, height=5.5)
par(mar=c(4.1,4.1,1.1,2.1))
pp.plotter.logit(model=results.mod1,x="Prejudice",Xc=Xc1,subtitle="Model 1",n=100000)
graphics.off()

### marginal effects / first differences ###
# create a function
fd.calculator <- function(model, Xc1, Xc2,seed=02143,n=10000){
  # require necessary libraries
  require(mvtnorm)
  # get necessary objects from our glm object
  coefs <- model$coefficients
#   varcov <- vcov(model)
  varcov <- weave(model)
  # get our simulated betas (beta tildes)
  set.seed(seed)
  sim.betas <- rmvnorm(n = n,
                       mean = coefs,
                       sigma = varcov)
  fd.ests <- (1 / (1 + exp(-Xc1 %*% t(sim.betas)))) - 
    (1 / (1 + exp(-Xc2 %*% t(sim.betas))))
  fd.ci <- c(quantile(fd.ests,.025),quantile(fd.ests,.9725))
  mean.fd <- mean(fd.ests)
  output.vec <- c(mean.fd,fd.ci)
  names(output.vec) <- c("First Diff", "Lower 95% CI", "Upper 95% CI")
  #return(list(output.vec,fd.ests))
  return(output.vec)
}
# create a function to add confidence bands
segmenter <- function(object,y,col="black"){
  segments(x0=object[2], x1=object[3], y0=y, col=col)
  segments(x0=object[2], y0=y-.08, y1=y+.08, col=col)
  segments(x0=object[3], y0=y-.08, y1=y+.08, col=col)
}

# calculate the weighted distribution of the prejudice variable for nh-white
# strong dems
mod1.svy <- svydesign(id=~1, weights=~weight, data=as.data.frame(model.1.X))
prej.prop.table <- as.data.frame(prop.table(svytable(~Prejudice,mod1.svy)))
for(i in 1:nrow(prej.prop.table)){
  prej.prop.table$Cum.Freq[i] <- sum(prej.prop.table$Freq[1:i])
}
prej.prop.table
#            Prejudice        Freq   Cum.Freq
# 1  0.333333333333333 0.015520012 0.01552001
# 2  0.416666666666667 0.037694840 0.05321485
# 3  0.458333333333333 0.052256378 0.10547123
# 4                0.5 0.428653694 0.53412492
# 5  0.541666666666667 0.116131156 0.65025608
# 6  0.583333333333333 0.107559393 0.75781547
# 7              0.625 0.087271234 0.84508671
# 8  0.666666666666667 0.047178833 0.89226554
# 9  0.708333333333333 0.027092719 0.91935826
# 10              0.75 0.014207547 0.93356580
# 11 0.791666666666667 0.045946745 0.97951255
# 12 0.833333333333333 0.007994807 0.98750736
# 13                 1 0.012492644 1.00000000

# from neutral (0.5) to 1
Xc.neutral <- c(1,apply(model.1.X[,2:13],2,mean))
Xc.neutral["Prejudice"] <- .5
Xc.1 <- c(1,apply(model.1.X[,2:13],2,mean))
Xc.1["Prejudice"] <- 1
fd.mod1.neutral.1 <- fd.calculator(model=results.mod1, Xc1=Xc.1, Xc2=Xc.neutral, n=100000)
# from neutral to .7083333 (around 92nd percentile)
Xc.point71 <- c(1,apply(model.1.X[,2:13],2,mean))
Xc.point71["Prejudice"] <-  .7083333
fd.mod1.neutral.point71 <- fd.calculator(model=results.mod1, Xc1=Xc.point71, Xc2=Xc.neutral, n=100000)
# from neutral to .625 (around 85th percentile)
Xc.point625 <- c(1,apply(model.1.X[,2:13],2,mean))
Xc.point625["Prejudice"] <-  .625
fd.mod1.neutral.point625 <- fd.calculator(model=results.mod1, Xc1=Xc.point625, Xc2=Xc.neutral, n=100000)
# plot these first differences
mod1.fd.ests <- c(fd.mod1.neutral.1[1], fd.mod1.neutral.point71[1], fd.mod1.neutral.point625[1])
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/mod1_marg_effects.pdf",
    width=7, height=3.5)
par(mar=c(4.1,4.1,1.1,2.1))
dotchart(x=mod1.fd.ests, xlim=c(-1,1), pch=19, cex=0.8,
         labels=c("Neutral to 1 (~p100)", "Neutral to .708 (~p89)", "Neutral to .625 (~p79)"),
         #main="Model 1: Marginal Effect of Prejudice",
         xlab="Marginal Effect")
segments(x0=0,y0=0,y1=4, lty=3, col="red")
segmenter(fd.mod1.neutral.1,1)
segmenter(fd.mod1.neutral.point71,2)
segmenter(fd.mod1.neutral.point625,3)
graphics.off()


###############
### model 5 ###
###############
### distribution of prejudice ###
# create a histogram of weighted and unweighted frequencies of prejudice
# run the function
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/mod5_distr_prejudice.pdf",
    width=7, height=4)
par(mar=c(4.1,4.1,1.1,2.1))
distr.prej.plotter(data=model.5.X,subtitle="Model 5")
graphics.off()

### predicted probabilities ###
# use the function to plot the predicted probabilities of turnout for values of prejudice
# set all variables to their mean
Xc.mod5 <- c(1,apply(subset(model.5.X[,2:15], model.5.X[,"Strong Dem"]==1),2,mean))
Xc.mod5["Strong Dem"] <- 1 # set strong dem =1 
# plot the predicted probability of prejudice
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/mod5_pred_probs.pdf",
    width=7, height=5.5)
par(mar=c(4.1,4.1,1.1,2.1))
pp.plotter.logit(model=results.mod5,x="Prejudice",Xc=Xc.mod5,interaction="Strong Dem x Prej",subtitle="Model 5",n=100000)
graphics.off()

### marginal effects / first differences ###
# from neutral (0.5) to 1
Xc5.neutral <- c(1,apply(subset(model.5.X[,2:15], model.5.X[,"Strong Dem"]==1),2,mean))
Xc5.neutral["Strong Dem"] <- 1
Xc5.neutral["Prejudice"] <- .5
Xc5.neutral["Strong Dem x Prej"] <- .5
Xc5.1 <- c(1,apply(subset(model.5.X[,2:15], model.5.X[,"Strong Dem"]==1),2,mean))
Xc5.1["Strong Dem"] <- 1
Xc5.1["Prejudice"] <- 1
Xc5.1["Strong Dem x Prej"] <- 1
fd.mod5.neutral.1 <- fd.calculator(model=results.mod5, Xc1=Xc5.1, Xc2=Xc5.neutral, n=100000)
# from neutral to .7083333 (around 92nd percentile)
Xc5.point71 <- c(1,apply(subset(model.5.X[,2:15], model.5.X[,"Strong Dem"]==1),2,mean))
Xc5.point71["Strong Dem"] <-  1
Xc5.point71["Prejudice"] <-  .7083333
Xc5.point71["Strong Dem x Prej"] <-  .7083333
fd.mod5.neutral.point71 <- fd.calculator(model=results.mod5, Xc1=Xc5.point71, Xc2=Xc5.neutral, n=100000)
# from neutral to .625 (around 85th percentile)
Xc5.point625 <- c(1,apply(subset(model.5.X[,2:15], model.5.X[,"Strong Dem"]==1),2,mean))
Xc5.point625["Strong Dem"] <-  1
Xc5.point625["Prejudice"] <-  .625
Xc5.point625["Strong Dem x Prej"] <-  .625
fd.mod5.neutral.point625 <- fd.calculator(model=results.mod5, Xc1=Xc5.point625, Xc2=Xc5.neutral, n=100000)
mod5.fd.ests <- c(fd.mod5.neutral.1[1], fd.mod5.neutral.point71[1], fd.mod5.neutral.point625[1])
# plot these first differences
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/mod5_marg_effects.pdf")
par(mar=c(4.1,4.1,1.1,2.1))
dotchart(x=mod5.fd.ests, xlim=c(-1,1), pch=19, cex=0.8,
         labels=c("Neutral to 1 (~p100)", "Neutral to .708 (~p89)", "Neutral to .625 (~p79)"),
         #main="Model 5: Marginal Effect of Prejudice",
         xlab="Marginal Effect")
segments(x0=0,y0=0,y1=4, lty=3, col="red")
segmenter(fd.mod5.neutral.1,1)
segmenter(fd.mod5.neutral.point71,2)
segmenter(fd.mod5.neutral.point625,3)
graphics.off()

###########################################
### models 1 & 5: plot marginal effects ###
###########################################
# create a plot
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/mod_1_5_marg_effects.pdf",
    width=7, height=4)
par(mar=c(4.1,4.1,1.1,2.1))
dotchart(mod1.fd.ests, xlim=c(-1,1), pch=19, cex=0.8,
         labels=c("Neutral to 1 \n (100th percentile)", 
                  "Neutral to .708 \n (92nd percentile)", 
                  "Neutral to .625 \n (85th percentile)"),
         xlab="Marginal Effect",
         #main="Marginal Effect of Prejudice", 
         lcolor="white")
# add a dashed line at 0
segments(x0=0,y0=0,y1=4, lty=3, lwd=2, col="gray48")
# plot the marginal effects from model 1
#points(x=mod1.fd.ests, y=1:3, pch=19, cex=.8, col="purple")
# plot the marginal effects from model 5
points(x=mod5.fd.ests, y=1:3-.25, pch=17, cex=.9, col="gray48")
# add a legend
legend(x="topleft",c("Model 1", "Model 5"), pch=c(19,17), 
       col=c("black","gray48"), cex=.8)
# plot the confidence intervals
segmenter(fd.mod1.neutral.point625,3,col="black")
segmenter(fd.mod5.neutral.point625,3-.25,col="gray48")
segmenter(fd.mod1.neutral.point71,2,col="black")
segmenter(fd.mod5.neutral.point71,2-.25,col="gray48")
segmenter(fd.mod1.neutral.1,1,col="black")
segmenter(fd.mod5.neutral.1,1-.25,col="gray48")
# add labels for at each marginal effect
text(x=mod1.fd.ests, y=1:3, labels=round(mod1.fd.ests,2), adj=c(1,-1), cex=.8, col="black")
text(x=mod5.fd.ests, y=1:3-.25, labels=round(mod5.fd.ests,2), adj=c(1,1.75), cex=.8, col="gray48")
graphics.off()



###############
### model 6 ###
###############
### distribution of prejudice ###
# create a histogram of weighted and unweighted frequencies of prejudice
# run the function
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/plots/mod6_distr_prejudice.pdf")
distr.prej.plotter(data=model.6.X,subtitle="Model 6")
graphics.off()

### predicted probabilities ###
# use the function to plot the predicted probabilities of turnout for values of prejudice
# set all variables to their mean
Xc.mod6 <- c(1,apply(model.6.X[,2:4],2,mean))
Xc.mod6["Strong Dem"] <- 1 # set strong dem =1
# plot the predicted probability of prejudice
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/plots/mod6_pred_probs.pdf")
pp.plotter.logit(model=results.mod6,x="Prejudice",Xc=Xc.mod6,interaction="Strong Dem x Prej",subtitle="Model 6",n=100000)
graphics.off()

### marginal effects / first differences ###
# from neutral (0.5) to 1
Xc6.neutral <- c(1,apply(model.6.X[,2:4],2,mean))
Xc6.neutral["Strong Dem"] <- 1
Xc6.neutral["Prejudice"] <- .5
Xc6.neutral["Strong Dem x Prej"] <- .5
Xc6.1 <- c(1,apply(model.6.X[,2:4],2,mean))
Xc6.1["Strong Dem"] <- 1
Xc6.1["Prejudice"] <- 1
Xc6.1["Strong Dem x Prej"] <- 1
fd.mod6.neutral.1 <- fd.calculator(model=results.mod6, Xc1=Xc6.neutral, Xc2=Xc6.1, n=100000)
# from neutral to .7083333 (around 92nd percentile)
Xc6.point71 <- c(1,apply(model.6.X[,2:4],2,mean))
Xc6.point71["Strong Dem"] <-  1
Xc6.point71["Prejudice"] <-  .7083333
Xc6.point71["Strong Dem x Prej"] <-  .7083333
fd.mod6.neutral.point71 <- fd.calculator(model=results.mod6, Xc1=Xc6.neutral, Xc2=Xc6.point71, n=100000)
# from neutral to .625 (around 85th percentile)
Xc6.point625 <- c(1,apply(model.6.X[,2:4],2,mean))
Xc6.point625["Strong Dem"] <-  1
Xc6.point625["Prejudice"] <-  .625
Xc6.point625["Strong Dem x Prej"] <-  .625
fd.mod6.neutral.point625 <- fd.calculator(model=results.mod6, Xc1=Xc6.neutral, Xc2=Xc6.point625, n=100000)
# plot these first differences
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/plots/mod6_marg_effects.pdf")
mod6.fd.ests <- c(fd.mod6.neutral.1[1], fd.mod6.neutral.point71[1], fd.mod6.neutral.point625[1])
dotchart(x=mod6.fd.ests, xlim=c(-1,1), pch=19, cex=0.8,
         labels=c("Neutral to 1 (~p100)", "Neutral to .708 (~p89)", "Neutral to .625 (~p79)"),
         xlab="Marginal Effect", main="Model 5: Marginal Effect of Prejudice")
segments(x0=0,y0=0,y1=4, lty=3, col="red")
segmenter(fd.mod6.neutral.1,1)
segmenter(fd.mod6.neutral.point71,2)
segmenter(fd.mod6.neutral.point625,3)
graphics.off()


#############################################
### simulate counterfactual turnout rates ###
#############################################
turnout.simulator <- function(model, seed=02143, n=10000, model5=0){
  # draw n sets of coefficients from the multivariate normal distribution
  coefs <- model$coefficients
  varcov <- weave(model)
  set.seed(seed)
  sim.betas <- rmvnorm(n = n,
                       mean = coefs,
                       sigma = varcov)
  # subset option to white, strong dems and add column of 1s
  if(model5==1){
    temp.data <- model$model[,2]
    wt <- model$model[,3][temp.data[,"Strong Dem"]==1]
    temp.data <- temp.data[temp.data[,"Strong Dem"]==1,]
    temp.data <- as.matrix(cbind(Intercept=1,temp.data))
    # create a new matrix for the no-prejudice counterfactual
    # that is, set all values of prej > 0.5 equal to 0.5
    temp.data.no.prej <- temp.data
    temp.data.no.prej[temp.data.no.prej[,"Prejudice"] > 0.5,"Strong Dem x Prej"] <- 0.5
    temp.data.no.prej[temp.data.no.prej[,"Prejudice"] > 0.5,"Prejudice"] <- 0.5
  }
  # take our data matrix and add a column of 1s for the intercept
  if(model5==0){
    temp.data <- as.matrix(cbind(Intercept=1,model$model[,2]))
    # create a new matrix for the no-prejudice counterfactual
    # that is, set all values of prej > 0.5 equal to 0.5
    temp.data.no.prej <- temp.data
    temp.data.no.prej[temp.data.no.prej[,"Prejudice"] > 0.5,"Prejudice"] <- 0.5
    wt <- model$model[,3]
  }
  # probability of turnout based on the model
  prob.turnout.prej <- (1 / (1 + exp(-temp.data %*% t(sim.betas))))
  prob.turnout.no.prej <- (1 / (1 + exp(-temp.data.no.prej %*% t(sim.betas))))
  # create a matrix to store predicted values
  pred.vals.prej <- apply(prob.turnout.prej, c(1,2), function(x){rbinom(n=1,size=1,prob=x)})
  pred.vals.no.prej <- apply(prob.turnout.no.prej, c(1,2), function(x){rbinom(n=1,size=1,prob=x)})
  # calculate distribution of turnout rates
  distr.turnout.prej <- apply(pred.vals.prej, 2, mean)
  distr.turnout.no.prej <- apply(pred.vals.no.prej, 2, mean)
  # calculate the distribution of turnout rates accounting for sampling weights
  wt.distr.turnout.prej <- apply(pred.vals.prej, 2, function(x) sum(x*wt)/sum(wt))
  wt.distr.turnout.no.prej <- apply(pred.vals.no.prej, 2, function(x) sum(x*wt)/sum(wt))
  # return the simulated distributions of turnout rates in a matrix
  wt.turnout.dist <- cbind(wt.distr.turnout.prej,wt.distr.turnout.no.prej)
  colnames(wt.turnout.dist) <- c("wt.turnouts.prej", "wt.turnouts.no.prej")
  return(as.data.frame(wt.turnout.dist))
}

# create a function to add confidence bands
segmenter2 <- function(object,y,col="black"){
  segments(x0=object[1], x1=object[2], y0=y, col=col)
  segments(x0=object[1], y0=y-.08, y1=y+.08, col=col)
  segments(x0=object[2], y0=y-.08, y1=y+.08, col=col)
}

# calculate the estimated turnout rate based on the data used in model 1
mean.data.turnout.prej <- svymean(~Turnout, mod1.svy)[1]
se.data.turnout.prej <- SE(svymean(~Turnout, mod1.svy))
t.crit <- qt(1-.025, df=degf(mod1.svy))
ci.data.turnout.prej <- c(mean.data.turnout.prej-se.data.turnout.prej*t.crit,
                          mean.data.turnout.prej+se.data.turnout.prej*t.crit)

# run the turnout simulator for model 1
mod1.turnout.dist <- turnout.simulator(model=results.mod1, seed=02143, n=10000)
# calculate the means and 95% CI
mean.mod1.turnout.prej <- mean(mod1.turnout.dist$wt.turnouts.prej)
ci.mod1.turnout.prej <- quantile(mod1.turnout.dist$wt.turnouts.prej,c(.025,.975))
mean.mod1.turnout.no.prej <- mean(mod1.turnout.dist$wt.turnouts.no.prej)
ci.mod1.turnout.no.prej <- quantile(mod1.turnout.dist$wt.turnouts.no.prej,c(.025,.975))

# run the turnout simulator for model 5
mod5.turnout.dist <- turnout.simulator(model=results.mod5, seed=02143, n=10000, model5=1)
# calculate the means and 95% CI
mean.mod5.turnout.prej <- mean(mod5.turnout.dist$wt.turnouts.prej)
ci.mod5.turnout.prej <- quantile(mod5.turnout.dist$wt.turnouts.prej,c(.025,.975))
mean.mod5.turnout.no.prej <- mean(mod5.turnout.dist$wt.turnouts.no.prej)
ci.mod5.turnout.no.prej <- quantile(mod5.turnout.dist$wt.turnouts.no.prej,c(.025,.975))

# plot the turnout rates and confidence intervals
means.prej <- c(mean.data.turnout.prej,mean.mod1.turnout.prej,mean.mod5.turnout.prej)
means.no.prej <- c(mean.mod1.turnout.no.prej,mean.mod5.turnout.no.prej)
# plot turnout means with prejudice
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/mod1_5_counterfactual_turnouts.pdf",
    width=7, height=4)
par(mar=c(4.1,4.1,1.1,2.1))
dotchart(x=means.prej, xlim=c(0,1), lcolor="white", xlab="Turnout Rate",
         pch=19, cex=0.8,labels=c("Estimated Turnout \n from the Data",
                                  "Predicted Turnout \n from Model 1",
                                  "Predicted Turnout \n from Model 5"))
# plot the turnout means with NO prejudice
points(x=means.no.prej,
       y=2:3-.25, pch=17, cex=.9, col="gray48")
# plot the confidence intervals
segmenter2(ci.data.turnout.prej,1)
segmenter2(ci.mod1.turnout.prej,2)
segmenter2(ci.mod5.turnout.prej,3)
segmenter2(ci.mod1.turnout.no.prej,2-.25,col="gray48")
segmenter2(ci.mod5.turnout.no.prej,3-.25,col="gray48")
# add a legend
legend(x="topleft",c("With Observed Prejudice", "Counterfactual with No Prejudice"), 
       pch=c(19,17), col=c("black","gray48"), cex=.8)
# add labels for at mean
text(x=means.prej, y=1:3, labels=round(means.prej,2), adj=c(.5,-1), cex=.8, col="black")
text(x=means.no.prej, y=2:3-.25, labels=round(means.no.prej,2), adj=c(.5,-1), cex=.8, col="gray48")
graphics.off()

##################################################################
### calculate the percentage of white, strong dems in the ANES ###
##################################################################
anes.svy <- svydesign(id=~1, weights=~weight, data=as.data.frame(anes))
pct.white.strong.dem <- as.data.frame(prop.table(svytable(~strong.dem.nh.white,anes.svy)))
#   strong.dem.nh.white      Freq
# 1                   0 0.9069892
# 2                   1 0.0930108



###################
### SENSITIVITY ###
###################
### excluding strong dems with prejudice<0.5 ###
# model 1: subset to strong dems (and whites)
model.1.X.nolt.5 <- na.omit(subset(model.1234.X,anes$strong.dem==1 & anes$nh.white==1 &
                                      anes$prejudice>=.5))
# model 2: subset to weak/lean dems (and whites)
model.2.X.nolt.5 <- na.omit(subset(model.1234.X,(anes$weak.dem==1 | anes$lean.dem==1) &
                              anes$nh.white==1 & anes$prejudice>=.5))
# model 3: subset to independents (and whites)
model.3.X.nolt.5 <- na.omit(subset(model.1234.X,anes$ind==1 & anes$nh.white==1 & 
                                      anes$prejudice>=.5))
# model 4: subset to republicans (and whites)
model.4.X.nolt.5 <- na.omit(subset(model.1234.X,(anes$strong.rep==1 | anes$weak.rep==1 |
                                            anes$lean.rep==1) & anes$nh.white==1 & anes$prejudice>=.5))
# model 5: all whites
model.5.X.nolt.5 <- na.omit(subset(model.5.X1,anes$nh.white==1 & anes$prejudice>=.5 ))

# model 6: all whites
model.6.X.nolt.5 <- na.omit(subset(model.6.X1,anes$nh.white==1 & anes$prejudice>=.5 ))

# run the function for each model
results.mod1.nolt.5 <- logit.Krup.Pist(model.1.X.nolt.5)
results.mod2.nolt.5 <- logit.Krup.Pist(model.2.X.nolt.5)
results.mod3.nolt.5 <- logit.Krup.Pist(model.3.X.nolt.5)
results.mod4.nolt.5 <- logit.Krup.Pist(model.4.X.nolt.5)
results.mod5.nolt.5 <- logit.Krup.Pist(model.5.X.nolt.5)
results.mod6.nolt.5 <- logit.Krup.Pist(model.6.X.nolt.5)

# output results to latex
results.list.nolt.5 <- list(results.mod1.nolt.5,results.mod2.nolt.5,results.mod3.nolt.5,
                              results.mod4.nolt.5,results.mod5.nolt.5,results.mod6.nolt.5)
se.list.nolt.5 <- lapply(results.list.nolt.5, weave.se)

# output clean results to LaTeX
stargazer(results.list.nolt.5,se=se.list.nolt.5,
          title=c("Sensitivity Test 2: Prejudice and White turnout in the 2008 election, by partisan group (ANES data)"),
          covariate.labels=cov.labels,dep.var.labels=c("Turnout"),column.sep.width="2pt")

### set prejudice =.5 for strong dems with prejudice<0.5 ###
# model 1: subset to strong dems (and whites)
model.1.X.set.5 <- na.omit(subset(model.1234.X,anes$strong.dem==1 & anes$nh.white==1))
model.1.X.set.5[model.1.X.set.5[,"Prejudice"] < 0.5,"Prejudice"] <- 0.5

# model 2: subset to weak/lean dems (and whites)
model.2.X.set.5 <- na.omit(subset(model.1234.X,(anes$weak.dem==1 | anes$lean.dem==1) &
                                    anes$nh.white==1))
# model 3: subset to independents (and whites)
model.3.X.set.5 <- na.omit(subset(model.1234.X,anes$ind==1 & anes$nh.white==1 & 
                                     anes$prejudice>=.5))
# model 4: subset to republicans (and whites)
model.4.X.set.5 <- na.omit(subset(model.1234.X,(anes$strong.rep==1 | anes$weak.rep==1 |
                                                   anes$lean.rep==1) & anes$nh.white==1 & anes$prejudice>=.5))
# model 5: all whites
model.5.X.set.5 <- na.omit(subset(model.5.X1,anes$nh.white==1))
model.5.X.set.5[model.5.X.set.5[,"Prejudice"] < 0.5 & model.5.X.set.5[,"Strong Dem"]==1,"Strong Dem x Prej"] <- 0.5
model.5.X.set.5[model.5.X.set.5[,"Prejudice"] < 0.5,"Prejudice"] <- 0.5

# model 6: all whites
model.6.X.set.5 <- na.omit(subset(model.6.X1,anes$nh.white==1))
model.6.X.set.5[model.6.X.set.5[,"Prejudice"] < 0.5 & model.6.X.set.5[,"Strong Dem"]==1,"Strong Dem x Prej"] <- 0.5
model.6.X.set.5[model.6.X.set.5[,"Prejudice"] < 0.5,"Prejudice"] <- 0.5

# run the function for each model
results.mod1.set.5 <- logit.Krup.Pist(model.1.X.set.5)
results.mod2.set.5 <- logit.Krup.Pist(model.2.X.set.5)
results.mod3.set.5 <- logit.Krup.Pist(model.3.X.set.5)
results.mod4.set.5 <- logit.Krup.Pist(model.4.X.set.5)
results.mod5.set.5 <- logit.Krup.Pist(model.5.X.set.5)
results.mod6.set.5 <- logit.Krup.Pist(model.6.X.set.5)

# output results to latex
results.list.set.5 <- list(results.mod1.set.5,results.mod2.set.5,results.mod3.set.5,
                            results.mod4.set.5,results.mod5.set.5,results.mod6.set.5)
se.list.set.5 <- lapply(results.list.set.5, weave.se)

# output clean results to LaTeX
stargazer(results.list.set.5,se=se.list.set.5,
          title=c("Sensitivity Test 3: Prejudice and White turnout in the 2008 election, by partisan group (ANES data)"),
          covariate.labels=cov.labels,dep.var.labels=c("Turnout"),column.sep.width="2pt")

### excluding strong dems with prejudice==1 ###
# model 1: subset to strong dems (and whites)
model.1.X.noprej1 <- na.omit(subset(model.1234.X,anes$strong.dem==1 & anes$nh.white==1 &
                                      anes$prejudice!=1))
# model 2: subset to weak/lean dems (and whites)
model.2.X.noprej1 <- na.omit(subset(model.1234.X,(anes$weak.dem==1 | anes$lean.dem==1) &
                                      anes$nh.white==1 & anes$prejudice!=1))
# model 3: subset to independents (and whites)
model.3.X.noprej1 <- na.omit(subset(model.1234.X,anes$ind==1 & anes$nh.white==1 & 
                                      anes$prejudice!=1))
# model 4: subset to republicans (and whites)
model.4.X.noprej1 <- na.omit(subset(model.1234.X,(anes$strong.rep==1 | anes$weak.rep==1 |
                                                    anes$lean.rep==1) & anes$nh.white==1 & anes$prejudice!=1))
# model 5: all whites
model.5.X.noprej1 <- na.omit(subset(model.5.X1,anes$nh.white==1 & anes$prejudice!=1))
# model 6: all whites
model.6.X.noprej1 <- na.omit(subset(model.6.X1,anes$nh.white==1 & anes$prejudice!=1))

# run the function for each model
results.mod1.noprej1 <- logit.Krup.Pist(model.1.X.noprej1)
results.mod2.noprej1 <- logit.Krup.Pist(model.2.X.noprej1)
results.mod3.noprej1 <- logit.Krup.Pist(model.3.X.noprej1)
results.mod4.noprej1 <- logit.Krup.Pist(model.4.X.noprej1)
results.mod5.noprej1 <- logit.Krup.Pist(model.5.X.noprej1)
results.mod6.noprej1 <- logit.Krup.Pist(model.6.X.noprej1)

# output results to latex
results.list.noprej1 <- list(results.mod1.noprej1,results.mod2.noprej1,results.mod3.noprej1,
                             results.mod4.noprej1,results.mod5.noprej1,results.mod6.noprej1)
se.list.noprej1 <- lapply(results.list.noprej1, weave.se)

# output clean results to LaTeX
stargazer(results.list.noprej1,se=se.list.noprej1,
          title=c("Sensitivity Test 1: Prejudice and White turnout in the 2008 election, by partisan group (ANES data)"),
          covariate.labels=cov.labels,dep.var.labels=c("Turnout"),column.sep.width="2pt")

### plot marginal effects for sensitivity tests
# sensitivity test 1
fd.sens1.mod1.neutral.71 <- fd.calculator(model=results.mod1.noprej1, Xc1=Xc.point71, Xc2=Xc.neutral, n=100000)
fd.sens1.mod5.neutral.71 <- fd.calculator(model=results.mod5.noprej1, Xc1=Xc5.point71, Xc2=Xc5.neutral, n=100000)
# sensitivity test 2
fd.sens2.mod1.neutral.71 <- fd.calculator(model=results.mod1.nolt.5, Xc1=Xc.point71, Xc2=Xc.neutral, n=100000)
fd.sens2.mod5.neutral.71 <- fd.calculator(model=results.mod5.nolt.5, Xc1=Xc5.point71, Xc2=Xc5.neutral, n=100000)
# sensitivity test 3
fd.sens3.mod1.neutral.71 <- fd.calculator(model=results.mod1.set.5, Xc1=Xc.point71, Xc2=Xc.neutral, n=100000)
fd.sens3.mod5.neutral.71 <- fd.calculator(model=results.mod5.set.5, Xc1=Xc5.point71, Xc2=Xc5.neutral, n=100000)
# plot these first differences
sens.mod1.fd.ests <- c(fd.sens1.mod1.neutral.71[1], fd.sens2.mod1.neutral.71[1], fd.sens3.mod1.neutral.71[1])
sens.mod5.fd.ests <- c(fd.sens1.mod5.neutral.71[1], fd.sens2.mod5.neutral.71[1], fd.sens3.mod5.neutral.71[1])
pdf("/Users/Daniel/Dropbox/gov2001_replication/shared/paper/tex_5-5-14/figures/sens_marg_effects.pdf",
    width=7, height=4)
par(mar=c(4.1,4.1,1.1,2.1))
dotchart(sens.mod1.fd.ests, xlim=c(-1,1), pch=19, cex=0.8,
         labels=c("Remove Prejudice =1", 
                  "Remove Prejudice < 0.5", 
                  "Set Prejudice < 0.5 to =0.5"),
         xlab="Marginal Effect",
         lcolor="white")
# add a dashed line at 0
segments(x0=0,y0=0,y1=4, lty=3, lwd=2, col="gray48")
# plot the marginal effects from model 5
points(x=sens.mod5.fd.ests, y=1:3-.25, pch=17, cex=.9, col="gray48")
# add a legend
legend(x="topleft",c("Model 1", "Model 5"), pch=c(19,17), 
       col=c("black","gray48"), cex=.8)
# plot the confidence intervals
segmenter(fd.sens3.mod1.neutral.71,3,col="black")
segmenter(fd.sens3.mod5.neutral.71,3-.25,col="gray48")
segmenter(fd.sens2.mod1.neutral.71,2,col="black")
segmenter(fd.sens2.mod5.neutral.71,2-.25,col="gray48")
segmenter(fd.sens1.mod1.neutral.71,1,col="black")
segmenter(fd.sens1.mod5.neutral.71,1-.25,col="gray48")
# add labels for at each marginal effect
text(x=sens.mod1.fd.ests, y=1:3, labels=round(sens.mod1.fd.ests,2), adj=c(1,-1), cex=.8, col="black")
text(x=sens.mod5.fd.ests, y=1:3-.25, labels=round(sens.mod5.fd.ests,2), adj=c(1,1.75), cex=.8, col="gray48")
graphics.off()



### excluding strong dems with prejudice==1

#na.omit(anes$V080001[anes$prejudice==1 & anes$strong.dem==1])
#anes.out <- anes[-1487,]
#anes.out <- anes.out[-1077,]
#na.omit(anes.out$V080001[anes.out$prejudice==1 & anes.out$strong.dem==1])


### MODELS 1-6 in TABLE 1 ###
# variables for models 1,2,3,4
#model.1234.X.out <- as.matrix(with(anes.out,cbind(
#  turnout,prejudice,age.0.1,educ.0.1,female,income.0.1,married,residence.0.1,
#  south,external.efficacy,interest,contact,close.race,weight,psu)))
#colnames(model.1234.X.out) <- c("Turnout","Prejudice","Age","Education","Female",
#                            "Income","Married","Years in residence","South",
#                            "Efficacy","Interest","Contacted by party",
#                            "Perceive election to be close","weight","psu")
# variables for model 5
#model.5.X.out <- as.matrix(with(anes.out,cbind(
#  turnout,prejudice,strong.dem,prejudice*strong.dem,age.0.1,educ.0.1,female,
#  income.0.1,married,residence.0.1,south,external.efficacy,interest,
#  contact,close.race,weight,psu)))
#colnames(model.5.X.out) <- c("Turnout","Prejudice","Strong Dem",
#                         "Strong Dem x Prej","Age","Education","Female",
#                         "Income","Married","Years in residence","South",
#                         "Efficacy","Interest","Contacted by party",
#                         "Perceive election to be close","weight","psu")
# variables for model 6
#model.6.X.out <- as.matrix(with(anes.out,cbind(
#  turnout,prejudice,strong.dem,prejudice*strong.dem,weight,psu)))
#colnames(model.6.X.out) <- c("Turnout","Prejudice","Strong Dem",
#                         "Strong Dem x Prej","weight","psu")



# model 1: subset to strong dems (and whites)
#model.1.X.out <- na.omit(subset(model.1234.X.out,anes.out$strong.dem==1 & anes.out$nh.white==1))
# model 2: subset to weak/lean dems (and whites)
#model.2.X.out <- na.omit(subset(model.1234.X.out,(anes.out$weak.dem==1 | anes.out$lean.dem==1) &
anes.out$nh.white==1))
# model 3: subset to independents (and whites)
#model.3.X.out <- na.omit(subset(model.1234.X.out,anes.out$ind==1 & anes.out$nh.white==1))
# model 4: subset to republicans (and whites)
#model.4.X.out <- na.omit(subset(model.1234.X.out,(anes.out$strong.rep==1 | anes.out$weak.rep==1 |
anes.out$lean.rep==1) & anes.out$nh.white==1))
# model 5: all whites
#model.5.X.out <- na.omit(subset(model.5.X.out,anes.out$nh.white==1))
# model 6: all whites
#model.6.X.out <- na.omit(subset(model.6.X.out,anes.out$nh.white==1))


# run the function for each model
#results.mod1.out <- logit.Krup.Pist(model.1.X.out)
#results.mod2.out <- logit.Krup.Pist(model.2.X.out)
#results.mod3.out <- logit.Krup.Pist(model.3.X.out)
#results.mod4.out <- logit.Krup.Pist(model.4.X.out)
#results.mod5.out <- logit.Krup.Pist(model.5.X.out)
#results.mod6.out <- logit.Krup.Pist(model.6.X.out)

#results.list.out <- list(results.mod1.out,results.mod2.out,results.mod3.out,results.mod4.out,results.mod5.out,results.mod6.out)
#se.list.out <- lapply(results.list.out, weave.se)

# output clean results to LaTeX
#stargazer(results.list.out,se=se.list.out,
#          title=c("Prejudice and White turnout in the 2008 election, with removals, by partisan group (ANES data)"),
#         covariate.labels=cov.labels,dep.var.labels=c("Turnout"),column.sep.width="2pt")


### prejudice scale
### post-election survey
prej.varlist.post <- colnames(with(anes,cbind(V085174a,V085174b,V085175a,V085175b)))
for(i in 1:length(prej.varlist.post)){
  anes[!is.na(anes[,prej.varlist.post[i]]) < 1,prej.varlist.post[i]] <- NA
}
# create the prejudice scale
lazy.scale.post <- with(anes,(V085174b - V085174a + 6)/12)
describe(lazy.scale.post)
intelligent.scale.post <- with(anes, (V085175b - V085175a + 6)/12)
describe(intelligent.scale.post)
anes$postjudice <- (lazy.scale.post + intelligent.scale.post) / 2
describe(anes$postjudice)

# calculate correlation between pre- and post-election measure
cor(anes$postjudice,anes$prejudice, use="na.or.complete")


### analysis:

### MODELS 1-6 in TABLE 1 ###
# variables for models 1,2,3,4
model.1234.X.post <- as.matrix(with(anes,cbind(
  turnout,postjudice,age.0.1,educ.0.1,female,income.0.1,married,residence.0.1,
  south,external.efficacy,interest,contact,close.race,weight,psu)))
colnames(model.1234.X.post) <- c("Turnout","Prejudice","Age","Education","Female",
                                 "Income","Married","Years in residence","South",
                                 "Efficacy","Interest","Contacted by party",
                                 "Perceive election to be close","weight","psu")
# variables for model 5
model.5.X.post <- as.matrix(with(anes,cbind(
  turnout,postjudice,strong.dem,prejudice*strong.dem,age.0.1,educ.0.1,female,
  income.0.1,married,residence.0.1,south,external.efficacy,interest,
  contact,close.race,weight,psu)))
colnames(model.5.X.post) <- c("Turnout","Prejudice","Strong Dem",
                              "Strong Dem x Prej","Age","Education","Female",
                              "Income","Married","Years in residence","South",
                              "Efficacy","Interest","Contacted by party",
                              "Perceive election to be close","weight","psu")
# variables for model 6
model.6.X.post <- as.matrix(with(anes,cbind(
  turnout,postjudice,strong.dem,prejudice*strong.dem,weight,psu)))
colnames(model.6.X.post) <- c("Turnout","Prejudice","Strong Dem",
                              "Strong Dem x Prej","weight","psu")

# model 1: subset to strong dems (and whites)
model.1.X.post <- na.omit(subset(model.1234.X.post,anes$strong.dem==1 & anes$nh.white==1))
# model 2: subset to weak/lean dems (and whites)
model.2.X.post <- na.omit(subset(model.1234.X.post,(anes$weak.dem==1 | anes$lean.dem==1) &
                                   anes$nh.white==1))
# model 3: subset to independents (and whites)
model.3.X.post <- na.omit(subset(model.1234.X.post,anes$ind==1 & anes$nh.white==1))
# model 4: subset to republicans (and whites)
model.4.X.post <- na.omit(subset(model.1234.X.post,(anes$strong.rep==1 | anes$weak.rep==1 |
                                                      anes$lean.rep==1) & anes$nh.white==1))
# model 5: all whites
model.5.X.post <- na.omit(subset(model.5.X.post,anes$nh.white==1))
# model 6: all whites
model.6.X.post <- na.omit(subset(model.6.X.post,anes$nh.white==1))

# run the function for each model
results.mod1.post <- logit.Krup.Pist(model.1.X.post)
results.mod2.post <- logit.Krup.Pist(model.2.X.post)
results.mod3.post <- logit.Krup.Pist(model.3.X.post)
results.mod4.post <- logit.Krup.Pist(model.4.X.post)
results.mod5.post <- logit.Krup.Pist(model.5.X.post)
results.mod6.post <- logit.Krup.Pist(model.6.X.post)

results.list.post <- list(results.mod1.post,results.mod2.post,results.mod3.post,results.mod4.post,results.mod5.post,results.mod6.post)
se.list.post <- lapply(results.list.post, weave.se)

# output clean results to LaTeX
stargazer(results.list.post,se=se.list.post,
          title=c(" Prejudice (post-election measure) and White turnout in the 2008 election, by partisan group (ANES data)"),
          covariate.labels=cov.labels,dep.var.labels=c("Turnout"),column.sep.width="2pt")

### prejudice scale
### average between pre- and post-election survey

# create the prejudice scale
lazy.scale.avg <- (lazy.scale+lazy.scale.post)/2
describe(lazy.scale.avg)
intelligent.scale.avg <- (intelligent.scale+intelligent.scale.post)/2
describe(intelligent.scale.avg)
anes$avgjudice <- (lazy.scale.avg + intelligent.scale.avg) / 2
describe(anes$avgjudice)

### analysis:

### MODELS 1-6 in TABLE 1 ###
# variables for models 1,2,3,4
model.1234.X.avg <- as.matrix(with(anes,cbind(
  turnout,avgjudice,age.0.1,educ.0.1,female,income.0.1,married,residence.0.1,
  south,external.efficacy,interest,contact,close.race,weight,psu)))
colnames(model.1234.X.avg) <- c("Turnout","Prejudice","Age","Education","Female",
                                "Income","Married","Years in residence","South",
                                "Efficacy","Interest","Contacted by party",
                                "Perceive election to be close","weight","psu")
# variables for model 5
model.5.X.avg <- as.matrix(with(anes,cbind(
  turnout,avgjudice,strong.dem,prejudice*strong.dem,age.0.1,educ.0.1,female,
  income.0.1,married,residence.0.1,south,external.efficacy,interest,
  contact,close.race,weight,psu)))
colnames(model.5.X.avg) <- c("Turnout","Prejudice","Strong Dem",
                             "Strong Dem x Prej","Age","Education","Female",
                             "Income","Married","Years in residence","South",
                             "Efficacy","Interest","Contacted by party",
                             "Perceive election to be close","weight","psu")
# variables for model 6
model.6.X.avg <- as.matrix(with(anes,cbind(
  turnout,avgjudice,strong.dem,prejudice*strong.dem,weight,psu)))
colnames(model.6.X.avg) <- c("Turnout","Prejudice","Strong Dem",
                             "Strong Dem x Prej","weight","psu")

# model 1: subset to strong dems (and whites)
model.1.X.avg <- na.omit(subset(model.1234.X.avg,anes$strong.dem==1 & anes$nh.white==1))
# model 2: subset to weak/lean dems (and whites)
model.2.X.avg <- na.omit(subset(model.1234.X.avg,(anes$weak.dem==1 | anes$lean.dem==1) &
                                  anes$nh.white==1))
# model 3: subset to independents (and whites)
model.3.X.avg <- na.omit(subset(model.1234.X.avg,anes$ind==1 & anes$nh.white==1))
# model 4: subset to republicans (and whites)
model.4.X.avg <- na.omit(subset(model.1234.X.avg,(anes$strong.rep==1 | anes$weak.rep==1 |
                                                    anes$lean.rep==1) & anes$nh.white==1))
# model 5: all whites
model.5.X.avg <- na.omit(subset(model.5.X.avg,anes$nh.white==1))
# model 6: all whites
model.6.X.avg <- na.omit(subset(model.6.X.avg,anes$nh.white==1))

# run the function for each model
results.mod1.avg <- logit.Krup.Pist(model.1.X.avg)
results.mod2.avg <- logit.Krup.Pist(model.2.X.avg)
results.mod3.avg <- logit.Krup.Pist(model.3.X.avg)
results.mod4.avg <- logit.Krup.Pist(model.4.X.avg)
results.mod5.avg <- logit.Krup.Pist(model.5.X.avg)
results.mod6.avg <- logit.Krup.Pist(model.6.X.avg)

results.list.avg <- list(results.mod1.avg,results.mod2.avg,results.mod3.avg,results.mod4.avg,results.mod5.avg,results.mod6.avg)
se.list.avg <- lapply(results.list.avg, weave.se)

# output clean results to LaTeX
stargazer(results.list.avg,se=se.list.avg,
          title=c("Prejudice (pre- & post-election average) and White turnout in the 2008 election, by partisan group (ANES data)"),
          covariate.labels=cov.labels,dep.var.labels=c("Turnout"),column.sep.width="2pt")
